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Abstract. Stellar clusters with the same general physical properties (e.g., total mass, age, and star-formation mode) may 
have very different stellar mass spectra due to the incomplete sampling of the underlying mass function; such differences are 
especially relevant in the high-mass tail of the mass function due to the smaller absolute number of massive stars. Since the 
ionising spectra of star-forming regions are mainly produced by massive stars and their by-products, the dispersion in the 
number of massive stars across individual clusters also produces a dispersion in the properties of the corresponding ionising 
spectra. This implies that regions with the same physical properties may produce very different emission line spectra, and 
occupy different positions in emission-line diagnostic diagrams. In this paper, we lay the bases for the future analysis of 
this effect by evaluating the dispersion in the ionising fluxes of synthetic spectra computed with evolutionary models. As 
an important consequence of the explicit consideration of sampling effects, we found that the intensities of synthetic fluxes at 
different ionisation edges are strongly correlated, a fact suggesting that no additional dispersion will result from the inclusion 
of sampling effects in the analysis of diagnostic diagrams; this is true for Hii regions on all scales, those ionised by single 
massive stars through those ionised by super stellar clusters. This finding is especially relevant, in consideration of the fact that 
real Hii regions are found in a band sensibly narrower than predicted by standard methods. Additionally, we find convincing 
suggestions that the He ii line intensities are strongly affected by sampling, especially during the WR phase, and so cannot be 
used to constrain the evolutionary status of stellar clusters. We also establish the range of applicability of synthesis models 
set by the Lowest Luminosity Limit for the ionising flux, that is the lowest limit in cluster mass for which synthesis models 
can be applied to predict ionising spectra. This limit marks the boundary between the situations in which the ionising flux is 
better modeled with a single star as opposed to a star cluster; this boundary depends on the metallicity and age of the stellar 
population, ranging from 10' to more than 10* Mq in the case of a single burst event. As a consequence, synthesis models 
should not be used to try to account for the properties of clusters with smaller masses. 
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1. Introduction and motivation 

In recent years a great deal of effort has been put into modeling 
more realistic ionising fluxes of ma ssive stars (e.g., the wm- 
BASic code by IPauldrach et all 1200 iL and references therein). 
These new model atmospheres, as included in evolutionary 
synthesis codes, are expected to produce more reliable ionis- 
ing spectra, which in turn should allow a better interpretation 
of the emission line spectra of ionised nebulae, th rough the use 
of photoionisation codes (e.g.. lSmith et alJ2002l) . 

Yet the ionising spectra of stellar clusters are produced 
by the most massive stars, whose absolute numbers can be 
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Strongly affected by the incomplete sampling of the underlying 
Initial Mass Function (IMF). The effects of incomplete sam- 
pling on the observed properties of star-forming regions can 
be cuiTently evaluated in two main ways: either by means of 
Monte Carlo synthesis models, or by means of an appropriate 
statistical formalism applied to the results of analytical synthe- 
sis models. By analytical synthesis models, or analytical mod- 
els for short, we will in the following refer to those popula- 
tion synthesis models computed by codes that neglect the issue 
of the incomplete sampling, and assume instead that the mass 
spectrum of stellar clusters can be represented by a continuous, 
analytical function that coincides with the underlying IMF. On 
the other hand, Monte Carlo codes simulate stellar populations 
by the stochastic generation of stars, which is stopped when ei- 
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ther a given number of stars or a given cluster mass is reached; 
the mass probability function for each star is given by the IMF, 
in such a way that the actual mass spectrum approaches the 
IMF as the number of stars increases, whereas large deviations 
occur for small clusters. This is assumed to be a more realistic 
representation of real stellar clusters than the one provided by 
analytical codes. 

The objective of this work is to study quantitatively the im- 
pact on emission line spectra of the sampling effects in the IMF, 
and to evaluate its consequences in the determination of the 
physical properties of Hii regions and galaxies. This is not a 
trivial problem. A simple, perhaps naive, approach would go 
through the following steps: (i) building synthetic clusters with 
different numbers of stars with a stochastic sampling of the 
IMF, ensuring in this way that only integer numbers of stars 
are included; (ii) obtaining the corresponding Spectral Energy 
Distributions (SEDs); (Hi) using the SEDs as inputs to a pho- 
toionisation code; and (iv) drawing conclusions based on the 
overall results. However, since the number of free parameters 
is pretty large, it is far more useful to split the problem into sev- 
eral smaller steps, each of which addressing a particular source 
of uncertainty. A possible procedure adopting this strategy is 
the following: 

- A first step is to evaluate quantitatively the sampling ef- 
fects on the ionising spectrum. Although the computation 
of photoionisation models is not involved in this step, it is 
important to characterize the resulting ionising flux with a 
view on the physics of photoionisation. 

- A second step is to determine how the different emission 
lines scale with the size (mass/number of stars) of the stel- 
lar cluster For example, the intensity of hydrogen recombi- 
nation lines roughly scale with the number of ionising pho- 
tons above the threshold energy, and hence with the number 
of massive stars in the cluster (for a given cluster age and 
given properties of the emitting gas). However, a similar 
scaling relation -provided it exists- has not yet been estab- 
lished for the case of forbidden emission lines. 

- A third step is to investigate the influence of both the stellar 
and the nebular properties on the results of photoionisation 
models. A large grid of photoionisation models would be 
required to explore exhaustively the parameter space, both 
along stellar population dimensions (i.e., IMF slope and 
mass limits, ages, stellar metallicity, model atmospheres, 
etc.), and along gas properties dimensions (density profile, 
chemical abundances, radius of the nebulae, covering fac- 
tor, etc.). In this step we can partially rel y on extensive stud- 
ies performed by different groups (e.g. Dopita et al. 2000, 
and references therein). Feedback from observations must 
also be taken into account. 

In this series of papers, we address each of these problems 
in turn. This first paper presents a quantitative evaluation of 
sampling effects on the predicted ionising flux for a wide range 
of input parameters in synthesis models, allowing us to obtain 
some conditions necessary, but not sufficient, to explain the 
correlations of observational data in emission-line diagnostic 
diagrams. In addition, we determine a mass limit under which 



stellar clusters should be better modeled with single individual 
stars rather than with stellar clusters. In a second paper we will 
determine the scaling relations of the emission-line spectrum 
with the initial mass/number of stars in the cluster The third, 
and last, paper of this series will be devoted to the analysis of a 
set of Monte Carlo simulations of evolutionary synthesis mod- 
els linked with a photoionisation code, where some specific 
properties will be assumed for the emitting gas. Throughout 
all these papers we will discuss the relationships between sam- 
pling effects and the ionising continuum, along with the result- 
ing emission-line spectra. 

The importance of observational feedback in this analysis 
cannot be overstated. A long-standing puzzle is posed by the 
positions in emission-line diagnostic diagrams of observed H ii 
regions ionised both by single stars and by stellar clusters cover 
an area narrower than the one expected on the basis of exten- 
sive grids of photoionisation models. For example, Dopita et alJ 
( 2000) built a grid of zero-age stellar clusters varying the metal- 
licity and the ionisation parameter of the model nebulae, and 
found that the area covered on diagnostic diagrams by these 
models is much broader than the area covered by observational 
data, suggesting the existence of a hidden parameter that con- 
strains the possible positions of real Hii regions in diagnos- 
tic diagiwis/niis^roblem has also been thoroughly discussed 
bv lStasiriska & Izotov.(,2003.) : these authors, who also consider 
the evolution of the ionising flux due to the aging of the stellar 
population, effectively highlight the difficulties of reproducing 
simultaneously all the observational constraints by means of 
traditional photoionisation models. 

Since the essence of this problem is finding a way of re- 
ducing the predicted dispersion, one would naively expect that 
the problem would worsen when the IMF sampling effects are 
taken into account: that is, that the positions occupied in diag- 
nostic diagrams by model nebulae would show a larger disper- 
sion as a consequence of taking into account the dispersion in 
the ionising flux predicted by stochastic models. Instead, we 
show here that the intensities of the predicted flux at different 
energies are strongly correlated, a fact implying that no addi- 
tional dispersion will result from the incomplete sampling of 
the IMF 

The structure of this paper is the following: in Section |2 
we review briefly the evaluation of sampling effects and the 
definition of the Lowest Luminosity Limit. The application of 
these concepts to the case of ionising spectra is presented in 
Section|3l Our main results are explained in Section|3 and the 
conclusions are presented in Section|5] 

2. Consequences of sampling effects in synthesis 
models 

In this section we describe how sampling effects in the IMF 
are evaluated. Some of these results have either been presented 
in previous papers, or can be found in books on statistics, e.g. 
Kendall & Stuart ( 1997) : nevertheless, we have chosen to sum- 
marize them here rather than in an appendix, because they are 
essential for the forthcoming discussion. Readers that are al- 
ready familiar with these tools may jump directly to Section|3| 
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2.1. The Lowest Luminosity Limit 

Sampling effects produce a dispersion in the results of synthesis 
models, and should therefore be taken into account when the 
properties of real clusters are interpreted by means of synthesis 
models. The relevance of such effects depends on the size of 
the system studied. A starting point to traduce this statement 
into quantitative terms is the definition of a criterion, related 
to the cluster size, allowing one to establish whether sampling 
effects in a particular cluster are playing a fundamental role 
or not. This alternative in turn would condition the mode of 
application of synthesis models: clusters small according to this 
criterion should necessarily be modeled with full consideration 
of sampling effects, whereas larger clusters could be modeled 
by means of traditional analytical methods. 

A suitable criterion to classify systems in this way 
is provided by the Lowes t Luminosity Limit (or LLL: 
ICervino & Lu ridiana'2003a'b^. The LLL is defined as the lumi- 
nosity of the most luminous individual star compatible with the 
physical properties (age, metallicity, etc.) of the cluster, and the 
criterion to follow is that an observed cluster whose luminos- 
ity is below the LLL cannot be modeled by means of a method 
that neglects sampling effects. The logic underlying this defini- 
tion is that an analytical model computed to reproduce a cluster 
below the LLL necessarily includes fractions smaller than one 
of the relevant stars, and therefore lack physical meaning. Note 
that the restriction imposed by this criterion is quite loose, since 
it constrains only those situations in which strong statistical ef- 
fects are undoubtedly present; but important statistical effects 
may well be present even in clusters above the LLL. 

Since the average luminosity A of a cluster scales linearly 
with the total mass Ai, the LLL is naturally associated to a 
lower mass limit Ai"™, which is the mass of a stellar system 
with a completely sampled IMF and luminosity equal to the 
LLL. 

The numerical value of the LLL depends on the isochrones 
and the atmosphere models used, does not depend on the IMF, 
and is only weakly dependent on the star formation history; 
on the contrary, it depends strongly on the age and metallic- 
ity of the synthetic cluster. Finally, the LLL obviously depends 
on the energy band in which the luminosity is defined. For 
a given luminosity band A the LLL defines an intrinsic limit 
below which a computation of the statistical dispersion asso- 
ciated to the incomplete sampling of the IMF is unavoidable 
for a meaningful application of synthesis models. Above the 
LLL, it is possible to provide a rough estimation of the dis- 
persion and a guideline t o determine the relevance of sam- 
pling effects: for example, 'Cervino & Luridiana' (2003a) esti- 
mate that the dispersion in the results of synthesis models for 
clusters with initial masses about 10 x Al™'" is equal or larger 
than 10% in the optical and infrared bands. Below the LLL, 
a proper statistical formalism is required to obtain quantitative 
estimations of the relative dispersion in the predicted quanti- 
ties; this topic will be covered in the next section. However, 
note that below the LLL the predicted average values of ob- 
servables that do not scale linearly with Ai, like logarithmic 
or rational functions of luminosities (i.e., equivalent widths or 
colours) can be severely biased with respect to actual obser- 



vations (see lCerviflo & Valls-Gabaudl2003l) . so that even a so- 
phisticated statistical formalism can fail to produce meaningful 
predictions. 

2.2. Estimation of sampling effects 

As stated in Section [2 sampling effects in synthesis models 
can be directly estimated by Monte Carlo methods, or alter- 
natively by a statistical formalism applied to the results of 
analytical synthesis models. A formalism of this kind, based 
on the original formalism proposed by Buzzoni ( 1989), has 
been developed in rece nt ye ars by jC ervino et al. (2000, 200 ih; 
ICervifio et ail ll2002bl) : and ICerviflo & Valls-Gabaud L2003I) . 
The method can be briefly summarised as follows. Let us 
assume that w, is the number of stars within a given mass 
range, normalised to the total mass; this value is given by the 
IMF and the star formation history. Assuming a Simple Stellar 
Population model (SSP, i.e. a model with coeval star formation, 
also called Instantaneous Burst model), w, is given by integrat- 
ing the IMF over the corresponding mass range. Each wi is ap- 
proximately the mean value of a Poisson distribution, and so 
the variance crj of each w,- is equal to w,. This assumption is a 
good approximation if th e mass interval that defi nes the w, val- 
ues is narrow enough: see lCervino et al.l(|2002bl) for a complete 
discussion of this topic. 

Given an observable property a, of a given star /, its con- 
tribution to the average normalised integrated property fi^^'' ob- 
tained by the synthesis model is given by w,a,, with a variance 
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af. The total variance of the average observable fi'^ 
is the sum of all the variances. Hence the relative dispersion is 
simply 



2^,1/2 



ssp 



1 



(1) 



where the last term gives the definition of Na introduced by 
I 1 1 1 ^ J 

IBuzzoni ( 1989). Since w; scales with the initial mass Ai trans- 
formed into stars in the cluster, the effective number of stars 
contributing to observable A also scales with the initial mass. 
This implies that, if the Ai value of an observed cluster is 
known, its predicted average A value is jUA - ^i^'' x Ai, and 
the dispersion around ^a is given by Na - N^'^'' x Ai. 

In addition, given another property b/, with contribution to 
the integrated average property jjg'' given by w,/?,- and variance 
(o"^-*"' ' possible to obtain the linear correlation coef- 

ficient between the observables A and B: 



p(A,B) 



cov(A,B) cov(A,B)' 



o-ao-b 



ssn ssp 



(2) 



where cov(A, B) is the covariance. Although cov(A, B) scales 
with Ai, p{A, B) is independent of it. The correlation coeffi- 
cient can also be obtained by a linear regression of A vs. B: a 
value of p(A, B) - \ implies that the two quantities are linearly 
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Fig. 1. Photoionisation cross sections ay(X') for the ions con- 
sidered in this paper, as obtained from the routine phf it . f in 
the photoionisation code Cloudy (.Ferland.1996) . 



resulting quantity, which we call the "effective rate of ionising 
photons" is defined by the expression; 



r 



hv 



ay(X') dv . 



(5) 



This quantity is obviously not an observable, but it will help 
to perform first-order estimations of the sampling effects on 
the ionising flux, via Eq.|3 We have computed Qa,,(X') for H°, 
He", He+, N", N+, 0°, 0+, S°, S+, C", Ne°, 
Ne+, Ar°, Ar+, and Ar++, although in this work flie Ar, Ne, 
and C ions will not be discussed'. The corresponding Oy values 
have been obtained from the subroutine phf it . f used for thi s 
purpose in the photoionisation code Cloudy (lFerlandlll996l) . 
and are shown in Fig.[n It is important to recall that since the 
recombination lines of X' arise from the recombination of X'^', 
they depend on Qa^{X') (e.g., Qo, (H'^) is related to the H i lines). 
On the other hand, since forbidden lines of X' arise from X' 
levels, they depend on Qo, ') (e.g. 2^, (O^) is related to the 
[Om] lines, and not to the [On] lines). 



dependent. In terms of a hypothetical diagnostic diagram of A 
vs. B, p(A, B) ~ 1 implies that the data points follow a narrow 
sequence. 

3. Methodology 

The simplest way to characterize the ionising flux of a stellar 
population is to consider both its intensity (in number of pho- 
tons emitted) and its shape (in terms of "effective rate of ion- 
ising photons" defined below). The emission rate of photons 
capable of ionising a particular ion (2(H°), 2(He^), etc.) with 
ionisation threshold frequency vq is given by: 



Qix') = 




As for the shape of the ionising flux, it is important to note 
that not all ionising photons are equivalent since the photoion- 
isation cross section depends strongly on v (see Fig. 0. The 
effectiveness of photons in ionising a ion X' is measured by the 
balance between X' and X'^' at a given point in the nebula: 

N(X''^^)Ne £°i4jrJv/hv)ay{X')dv 
N(X') ~ ac(X', T) ' 

where A^(X'+'), N{T), and A^^ are the number densities of ions 
X\ and electrons respectively, ac{X', T) is the recombi- 
nation coeflicient of the ground level of X'^' to all levels of 
X', ay{X') is the absorption cross section of the ion X' from 
its grou nd level with thre shold v,, and AnJy is the local flux 
density llOsterbrocklfl989h . Since we are ultimately interested 
in determining the emission spectra, which in turn depend on 
the relative ionic abundances, it seems physically meaningful 
to weigh the emission rate of photons capable to ionise a par- 
ticular ion by the photoionisation cross section for that ion, in 
order to describe the luminosity actually seen by each ion; The 



4. Results 

In this section we will describe our results. All the computa- 
tions assume an instantaneous burst of star formation and a 
ISalpeter ( 1955) IMF in the mass range 0.1 to 120 Mq. 

To assess the possible systematic effects associated with 
the choice of stellar tracks, model atmospheres, interpolation 
schemes, etc, we have used two diff'erent evolutionary synthe- 
sis codes: 

- The code by ICervifio et aP (l2002al) (hereinafter CMHK), 
which adopts evolutionary tracks with standard mass- 
loss rates ( S challer et al. 1992; Schaerer et al. 1993^3 
Charbonnel et al.l tl993l) and the model atmospheres by 
Schaerer & de KoteJ ( Il997|) (CoStar) for main sequence 
hot st ars more massive th an 20 Mp, by [Schmutz e t alJ 
( 1 19921) for WR stars, and bv lKuruczl ( ll99l]) for the remain- 
ing stars. 

- A modified version of Starburst99 llT,eithereretal.ll999h . 
where the main modification is the computation of the 
QoyiX') values and the corresponding minimum masses. 
The use of Starburst99 is only for comparison purposes, 
since we are intereste d in studying the i mpact of the new 
atmosphere models bv lSmith et al J (12002 ) for WR and hot 
stars, implemented in this code. Note that Starburst99 as- 
sumes Blackbody spectra for stars with log T^fj > 4.778 or 
log Teff < 3.3 except for WR stars (defined by the stellar 
temperature, log T, > 4.4, the hydrogen surface abundance 
Xs < 0.4, and the minimum mass for WR fo rmation). The 
stellar tracks used have high mass-loss rates JMevnet et alJ 
1994). 

In both codes, the spectrum associated to a given star is 
given by the closest model atmosphere found either in the 
logg - Teff plane or in the g - Tejf- For the computation of 
th e isochrones, both co des follow the prescriptions described 
inlCerviiioetal](l200ll). 

' The complete set of results can be found at the CMHK web pages 
at .http://laeff ■ inta.es/users/mcs/SED/.. 
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4.1. Effective rate of ionising photons QaSX') 

The time evolution of Qa^,{X') predicted by Starburst99 and 
the CMHK code is shown in Fig. Hand Irrespectively, where 
the five main panels correspond to the metallicities indicated^. 
Each main panel is further subdivided into four sub-panels, 
each showing the results for a different set of ions, coded ac- 
cording to the lines and symbols indicated. The top-left panel 
corresponds to H", He", and H^; the top-right panel to N", N^, 
and the bottom-left panel to O", 0+, and 0++; and the 
bottom-right panel to S", S+, and 

There are striking differences between the Qa,{X') be- 
haviour predicted by the two codes, a fact that emphasises the 
importance of the models' input - evolutionary tracks and at- 
mosphere models - as a source of systematic eff'ects. An exam- 
ple is the difference in the bumps of the 2fly(He^) curves: these 
bumps reveal the life cycle of WR stars, which provide a sub- 
stantial amount of hard photons. The low Z models computed 
with Starburst99 show such bumps at about 3 Myr, whereas 
the CMHK models do not; the difference is mainly due to the 
different evolutionary tracks, because enhanced mass-loss rates 
produce larger numbers of WR stars. A further difference is the 
larger 2n,,(He^) value at high Z in the CMHK code as com- 
pared to Starburst99. This difference can be ascribed to differ- 
ences in the atmosphere models implemented in the two codes: 
in general, the CMHK code produces a harder ioni sing flux due 
to the use of the CoStar and .Schmutz et alJ (Il992h atmosphere 
models. This effect is larger for high metallicities during the 
WR phase. The impact of different atmosp here models in syn - 
thesis codes has been recently discussed bv lSmith et alJ(l2002h . 
and we refer to that paper for further details. 

There are two other interesting features in Fig.|3and|3that 
deserve discussion: 

1. Although the behaviours of ^^^(He"^) and Qa^O^^) are 
remarkably similar, as a consequence of their ionisation 
edges being almost identical (54.418 eV for He^ and 
54.936 eV for O^^), there are slight differences between 
them due to the greater sensitivity to low energy photons of 
av(He+) with respect to flv(0++). It follows that 2fl.(0+^) 
decreases more rapidly than 2av(He^) when the hard ionis- 
ing flux decreases, that is for ages f > 4 Myr. 

2. The tails of the (He^) and Qa,{0^^) curves during the 
post-WR phases show an interesting behaviour in both sets 
of results. For example, at Z=0.008 there is a 2a,,(He^) 
peak around 6 Myr in the case of stellar tracks with en- 
hanced mass-loss rates (Fig.|2}, and a peak around 6.5 Myr 
in the case of stellar tracks with standard mass-loss rates 
(Fig. |3- Yet, although the bumps in these curves reveal 
the presence of WR stars in the stellar population, the WR 
phase ends at 5.9 Myr in the case of enhanced mass-loss 
rates, and at 4.7 Myr in the case of standard mass-loss rates. 
The explanation of this apparent contradiction lies in the 
track interpolation technique used to obtain the isochrones: 
even though it is assumed that the evolutionary tracks fol- 
low a continuous sequence, there is a discontinuity in the 

^ Note the different scale along the 2a,,-axis for Z=0.020 and 
Z=0.040. 
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Fig. 4. Evolutionary tracks ( dashed Une) and isoc hrones (solid 
lines) from Starburst99 an d iMevnet et al]lll994 . 

stellar tracks between WR and non-WR phases. The situ- 
ation is illustrated in Figure^] where isochrones for 4.5, 6 
and 6.3 Myr obtained from Starburst99 are compared with 
the evolutionary tracks for 60, 40 and 25 Mq with enhanced 
mass-loss rates. Note that some stars in the 6 Myr isochrone 
(where no WR stars should be present) reach log T^fj larger 
than 4.4, and they are therefore (wrongly) assigned to the 
WR type. A similar problems occurs in the CMHK code, 
and in general in all synthesis codes, since continuity in the 
tracks is implicitly ass umed (see, however, an alter native 
approach presented in ICervifio & Mas-Hessel[l994l) . The 
isochrones are mathematically consistent, but this is not 
a guarantee that they are physically consistent. This prob- 
lem is amplified in the tracks with enhanced mass-loss rates 
since adjacent tracks are far more different than the corre- 
sponding tracks with standard mass-loss rat es. This may su - 
perficially be considered a technical detail llMaedeil2002l) . 
but since it may produce unphysical results, it is worth in- 
vestigating it in detail (Cervino 2003, in preparation). 



4.2. Minimum masses M"""[QaXX')] 

In analogy to our definition of the minimum mass associated to 
the LLL, we will define in the following a minimum mass asso- 
ciated to the effective rate of ionising photons, M"""[Qa-,iX')], 
which is defined as the mass of a stellar system with a com- 
pletely sampled IMF and total effective rate of ionising photons 
equal to the maximum of all the effective rates provided by the 
individual stars in the system. 

The evolution of M"™[Qay(X')] for the ions selected be- 
fore, computed with Starburst99 and the CMHK code, are 
shown in Fig.|5]and|S]respectively. The organisation of the fig- 
ures is the same as in Fig. |3 Perhaps the most striking fea- 
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ture of these results is the huge mass range covered, which 
spans nearly three orders of magnitude, from about 3 x 10^ 
to more than 10^ Mq for some of the important ions, indepen- 
dently of the code used. As expected, the larger the ionising po- 
tential, the larger the minimum mass M'"'"[Qa^iX')], since in- 
complete sampling affect mainly the most massive stars, which 
are those contribute most to this property. Likewise, S° has the 
lowest Ai"""[Qa^(X')] value: the ionising potential of S" is only 
10.36 eV, so the effective rate of ionising photons for this ion is 
scarcely affected by sampling effects. On the other hand, He^ 
is the ion more affected by sampling effects: this result con- 
firms that He ii lines cannot be used as a strong constrain on the 
properties of stellar clusters, as suggested by iLuridiana et alJ 
(2003b). 

The most prominent peak in Ai'"'"[Qa^,{X')] corresponds 
to the WR phase of the cluster In the case of Starburst99, 
there is a tendency to obtain larger M"""[Qa,.iX')] values 
during the WR-phase at intermediate metallicities (Z=0.008), 
and an asymmetric "U-shaped" curve at low metallicities. 
We have tested this behaviour obtaining the corresponding 
M'"'"\Qa (X')] with the CMH K code using the evolutionary 
tracks bv lMevnet et a and the "U-shaped" behaviour 

remains. However, M"""[Qa^(X')] tends to be lower for larger 
metallicities. For the case of the CMHK code, the effect is ex- 
plained easily: since the WR atmosphere models are the same 
for all the metallicities, the larger the metallicity, the larger 
the mass-loss rate, the lower the luminosity, and the lower 
M"'"[Qa,(X')]. In the case of Starburst99 the effect of the 
metallicity-dependent WR atmosphere models must also be 
taken into account together with the variation of the luminosity 
with metallicity. In fact, examining the ionising fluxes quoted 
for the WR stars in tables 3 and 4 of Smith et al. (2002), it is 
found that the flux of ionising photons is not a monotonic func- 
tion of the metallicity: the largest values for the rate of ionis- 
ing photons correspond to some of the WR models at Z=0.008 
(Z=0.4 Zq in their paper). This would be the origin of the larger 
At"""[2fl^(X')] value at this metallicity. This also shows the im- 
pact of the use of metallicity-dependent atmosphere models for 
WR stars. 

The "U-shaped" behaviour during the WR-phase is due to 
the evolutionary tracks used: since the lower the number of 
stars that dominate the luminosity, the larger the Al""" [Qa^X')] 
value, its value depends on the relative rate of ionising pho- 
tons emitted by WR stars with respect to the other hot star s 
in the synthetic cluster. In the tracks by iMevnetet alJ lll994l) . 
the WR stars at the beginning of the WR phase are also the 
most luminous stars in the cluster (see, for example, the 60 
track in Fig.|3 with the harder ionising flux: this fact produces 
the first peak in Ai"""[Qa^,(X')]. When the system evolves, the 
luminosity of WR stars decreases and the other non-WR, hot 
stars give only a small contribution (see, for example, the 40 
Mq ti-ack in Fig.EJ, leading to the decrease in M'"'"[QaJX')]. 
However, since the non-WR stars become cooler as the age in- 
creases, at the end of the WR-phase WR stars become again the 
almost only source of ionising flux, producing the secondary 
peak. Note also that, for Z=0.040, the end of the WR-phase is 

^ Which results are also available at the CMHK web pages. 



marked by the highest peak in the curve: this behaviour is due 
to the longer duration of the WR-phase at this metallicity and 
to the dependency of the lifetime of massive stars with mass, 
which in these tracks is not monotonic: stars with initial mass 
of 85 and 120 Mq stars end their evolution as WR at 5 and 8 
Myrs respectively. As a consequence, in such an age range and 
at this metallicity, the ionising flux is dominated by the less 
populated tail of the IMF, so that sampling effects dominate 
and a larger M"""[Qa,,(X')] value is obtained. 

The same general explanations also apply to the case of 
the CMHK results with standard mass-loss rate. However, for 
these tracks WR and non-WR stars have more similar lumi- 
nosities and the first peak does not appear (the curve is "J- 
shaped"). In this case the maximum in M"""[Qay(X')] occurs 
at Z=0.004, because the WR phase is very short and almost ab- 
sent at Z=0.001 . In fact, at this metallicity, there no star ends its 
evolution as a WR; the WR-phase only appears in the middle 
of the 120 Mq track, which ends as a supergiant. 

There are also several peaks in M™{QaSX')\ at the oldest 
ages considered in these plots. These features are an artefact of 
the way in which model atmospheres are assigned to stars in 
the synthetic HR diagram: since the smooth evolutionary path 
of each star in the HR diagram is coupled to a discrete ensemble 
of atmosphere models, a peak is produced when stars suddenly 
switch from one atmosphere model to the next. In the case of 
Starburst99, the effect is not too strong because this code uses 
a finer atmosphere grid. In the case of the CMHK code, the 
switch from the CoStar to the .Kurucz C1991 ) atmosphere mod- 
els for stars in the main sequence around M — 20 M© produces 
the additional bumps in M"""[Qa^.{X')] with the most promi- 
nent ones ending at t k 9.4 Myr for Z^O.OOl and Z=0.004, t ^ 
9 Myr for Z=0.008, f ^ 8.1 Myr for Z=0.020, and f * 7 Myr 
for Z=0.040, corresponding to the turn-off age of a 20 M© star 
(point 13 in the tracks), the phase at which CoStar models are 
no longer used. Additional small bumps are found at the turn- 
off age of a 25 M©. We want to remark again that these bumps 
are model artifacts, and that they would not appear if the atmo- 
sphere models were interpolated, instead of being computed by 
means of the "closest atmosphere model" approach. 

Note that M"""[Qa^{X')] corresponds to the amount of gas 
transformed into stars, and not to the total mass of the cluster 
(stars plus gas, that is dynamical mass, which would be approx- 
imately a factor of ten larger), and that a Salpeter IMF in the 
mass range 0. 1 to 120 M© has been assumed. Table^gives con- 
version factors for different mass limits and also the numbers 
of stars for a Salpeter IMF. The results of Starburst99 are usu- 
ally normalised to 10* Mq in the mass range 1-120 Mq, while 
CMHK normalise the results to 1 M© in the mass range 2-120 
Mo. 

Using the values in Table [T] a value of M"""[QaJX')] 
around 3 x 10^ implies that there are, on average, two stars 
in the mass range 40-120 Mq, or 21 stars in the mass range 
8-120 Mq. Since the most massive stars are also the most lu- 
minous ones, the ionisation of the cluster is accounted for by a 
few stars. The presence of a few WR stars, with a harder ion- 
ising flux, produces the increases in M"""[Qa,.{X')], especially 
for ions with higher energy edges. A straightforward conclu- 
sion follows from this result: the ionisation structure for clus- 
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IMF mass range 






Quantity 


0.1-120 


1-120 


2-120 8-120 


25-120 


40-120 


Mass into stars 
Number of stars 


1 

2.8290 


0.3962 
0.1262 


0.2912 0.1442 
0.0494 0.0074 


0.0668 
0.0014 


0.0428 
0.0007 



Table 1. Conversion factors for cluster masses and numbers of stars for different IMF mass ranges assuming a Salpeter slope. 



ters with Ai close to Ai"""[Qa^iX')] would be better reproduced 
by the ionising spectrum of a single star rather than by the syn- 
thetic spectrum produced by a synthesis code, which always 
includes an ensemble of stars. 

It is necessary to recall, in this context, that M"""[Qa^(X')] 
is, in fact, a rather restrictive lower limit for the use of inte- 
grated ionising spectra obtained by synthesis models. A better 
way to estimate when a composite spectrum begins to be a bet- 
ter approximation is given by the relative dispersion in Qa,.(X') 
corresponding to Al = Al"™[2a^(X')], which is shown in 
Fig-H As expected, the dispersion varies strongly with the age, 
metallicity and ion under consideration, but, in general, it has a 
value around 60% (or, equivalently, N ~ 3). This means that in 
order to obtain a dispersion smaller than 10% {N = 100), clus- 
ters with initial masses M > 33 x M'"'"[Qa,(X')] are needed. 
This is a good lower limit to ensure that the use of the compos- 
ite ionising spectrum obtained from synthesis models is appro- 
priate. For clusters with masses in the range [Ai""", 33 x Ai"""] 
the situation is fuzzier, and any result obtained from codes 
which use a fully-sampled IMF must be taken with extreme 
caution. The only proper alternative solution in this case is to 
use Monte Carlo simulations. 

4.3. Correlations between ionising bands 

In this section, we will discuss the information that can 
be obtained through the study of the correlation coefficients 
plQciyiX'), QciyiYO] between the effective rate of photons of the 
ions X' and Y^. For the sake of simplicity, we will use in the 
following the notation p(X', Y^) = p{Qa,.{X'), 2„„(F^)). 

As indicated in Section a linear correlation coefficient 
p{X', Y') close to unity indicates that Qa.iX') and Qa,{Y') are 
essentially produced by the same stars within a stellar popu- 
lation. Note that the correlation between observable quantities 
is not affected by sampling effects in the IMF; nevertheless, 
we take advantage of the statistical formalism introduced in 
Section IT2I for the analysis of the effects of incomplete sam- 
pling to compute the values of the correlation coefficients. 

The correlation coefficients between effective rates of ion- 
ising photons may provide useful insights into the concept of 
ionisation correction factors (icfs) used in abundance determi- 
nations. A correlation coefficient close to one is a necessary 
condition to relate the abundance of two different ions in a 
ionised nebula independently of the ionising stellar population. 
The assumption of the existence of general recipes relating the 
abundances of two ions is conceptually analogous to the one 
underlying the definition of the icfs. However, a correlation co- 
efficient close to one is not a sufficient condition for a correct 
definition of the icfs, since the properties of the gas might also 



be relevant in t he determinatio n of the ionisation structure of a 
nebula (Luridia na et al.l2003ah . 

The correlation coefficients between the effective rates of 
ionising photons of different ions may also be used to make pre- 
dictions on diagnostic diagrams. However, in the case of diag- 
nostic diagrams an additional difficulty appears, because they 
do not relate absolute intensities but ratios. Hence, e.g., the cor- 
relation coefficients p(H'',N'') and p(H",0^) are not sufficient 
to derive a correlation in the [Nii] 6584/Ha' vs. [Oiii] 5007/Hj8 
diagram, since it is also necessary to know p(N°, O^). In addi- 
tion, even in the case of a full linear correlation {p - \), the 
tilt of the regression line and the position of the mean value are 
required, so a photoionisation model is always needed. For the 
time being, we will only obtain the correlation coefficients, to 
see whether correlations exist. In Paper III we will apply these 
results to make predictions on the dispersion in the diagnostic 
diagram. 

The results of the calculation of the correlation coefficients 
are shown in Fig.|8lfor selected ions and diagnostic diagrams. 
It is interesting to compare this figure with Fig. |6l where 
M"'"{QaSX')'\ is given as a function of time for different metal- 
licities and ions. A similar behaviour is observed in both fig- 
ures: the larger the difference of the Ai'"'"[Qay{X')] values, the 
lower the correlation coefficient. This is easily understood: a 
larger difference in Ai""" means that the stars that contribute 
the most to the corresponding effective rates of ionising pho- 
tons are different, and hence, the effective rates are loosely cor- 
related. 

In general, all correlation coefficients are close to unity un- 
til the WR phase begins. This phase introduces extra ionising 
photons associated with a handful of luminous stars, and hence 
destroy the correlations between the effective rates of ionis- 
ing photons. Correlation coefficients closer to one are prefer- 
entially found at low metallicities, due to the paucity of WR 
stars. 

Fig.|8]is divided into five main panels corresponding to dif- 
ferent metallicities. Each main panel is further subdivided into 
four sub-panels: the top-left sub-panel shows p(H'',He''), and 
P(H'', He^), the top-right sub-panel shows the correlation coef- 
ficients of ions with similar ionising potential, the bottom-left 
sub-panel shows the correlation coefficients related to the 7/ pa- 
rameter, and the bottom-right sub-panel the correlation coeffi- 
cients related to the [N 11] 6584/Hq' vs. [O m] 5007/Hj6 diagram. 

The p(H°, He") and p(H'', He^) correlation coefficients have 
values significatively different from 1 starting from the on- 
set of the WR-phase. The larger the metallicity, the lower the 
correlation coefficient, implying that the ratios 2(He'')/2(H") 
and 2(He^)/2(H°) are not robust constraints of the gen- 
eral properties of a cluster (hke the upper mass limit in the 
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IMF, the IMF slope, and the age): they are quite dependent 
on the particular stellar populations in the cluster and their 
value can strongly vary from cluster to cluster, especially for 
poorly populated high metallicity regions (see, for example, 
[Sresolin & Kennicutt 2002). On the other hand, for the same 
reason they can be used to infer the relative abundance of par- 
ticular populations in a given observed cluster 

The correlation coefficients shown in the top-right sub- 
panels are pCH", 0°) (the two ions having ionising potentials 
of 13.599 and 13.618 eV respectively), p(He°,S+) (ionising 
potentials of 24.588 and 23.330 eV respectively), p(He+, 0++) 
(ionising potentials of 54.418 and 54.936 eV respectively), and 
p(0^,S^^) (ionising potentials of 35.118 and 34.830 respec- 
tively). With the only exception of p(He^, O^^), the correlation 
coefficients are pretty close to 1 in all the cases, as expected 
during the WR phase. In contrast, the correlation coefficient of 
p(He^,0^^) does not have a constant value, but surprisingly 
varies between values close to one and values close to zero. 
The zero value does not necessarily imply that the two quanti- 
ties are not correlated in the corresponding age ranges: it rather 
indicates that one or both of the go, values are zero, that is 
p(He^, O^^) = implies the absence of photons above 54.418 
eV in the model atmospheres of stars at the corresponding ages. 

The bottom-left sub-panels correspond to correlation coef- 
ficients rel ated to the y? parameter, that is the ratio of O^/O^^ 
to S+/S++ jVflchez & Pagellll98a) . Note that the 77 parameter 
is measured from forbidden lines, thus it is related to Qa,.{0^^), 
Gi/,(0^), Qfl.XS") and gfl„(S+). In general, the correlation co- 
efficients involving are markedly lower than 1 . Remember 
that S" was the ion with the lowest Al""", which also means 
that it is produced by a population of stars different from those 
responsible for the ionising flux (in fact, it is produced by the 
ionising stars plus other stars): this fact naturally produces a 
poor correlation. As in the case of p{H°, He°) and p(H°, He+), 
the T] parameter is not a good parameter to constrain the age of 
the cluster, but it is a good parameter to determine the relative 
proportion of non-ionising and ionising stars, and the softness 
of the ionising radiation for specific clusters. 

Finally, the interesting correlation coefficients p(H'^,0^) 
and p(N'',0^), together with p(H°,N°), can be used to esti- 
mate the dispersion in the [Nn] 6584/Hq' vs. [Oiii] 5007/H/3 
diagram. The effective rates of ionising photons Qa^.iH^) and 
Gi/,.(N") are strongly correlated, implying that the abundances 
of the two ions are linearly dependent. In turn, this relation sug- 
gests that both elements will show a similar correlation with 
O^, a fact confirmed by piVf', O^) presenting a behaviour sim- 
ilar to p(N'^, O^). This has an additional implication: the po- 
sition of data points in the [Nii] 6584/Hq' vs. [Oiii] 5007/H;8 
diagram depends on two variables only (either O^ and H" or 
O^ and N"), and so the data points will lie on a uniparamet- 
ric line rather than being spread on the entire plane. Hence, 
the correlation coefficients obtained through the analysis of the 
sampling effects may explain the fact that observational data 
follow a narrow sequence. However, the sequence is unipara- 
metric only for a given metallicity, and it could be different 
when metallicity varies, yielding a larger scatter in the diagnos- 
tic diagram. That is, in fact, the situation in the lower region of 



this diagram, which is populated by data points corresponding 
to different metallicities. 

5. Conclusions 

In this paper we have presented an evaluation of the dispersion 
in the properties of the ionising flux obtained from synthesis 
models due to the IMF sampling. The ionising fluxes have been 
characterised by appropriate physical quantities, the effective 
rates of ionising photons, which are obtained from the integra- 
tion of the ionising flux weighed by the photoionisation cross 
section of different ions, with the future goal of using them in 
the evaluation of sampling effects on emission line spectra. 

We have obtained the effective rate of ionising photons for a 
wide set of metallicities and atmosphere models, comparing the 
results of two different codes to assess some of the systematic 
effects. The analysis of the effective rates obtained with differ- 
ent model atmospheres is consistent with previous studies on 
the impact of the atmosphere models on the ionising flux of star 
forming regions, and it additionally provides a proper statisti- 
cal ground for such assessments. We have identified problems 
in the computation of the ionising flux in regimes where stel- 
lar evolution does not present a continuous behaviour with the 
initial stellar mass, such as at the end of the WR phase. 

We have computed the corresponding minimum initial 
cluster masses for which the ionising continuum obtained by 
synthesis models can be applied. The minimum masses range 
from 10^ to more than 10^ M© depending on the metallicity 
and the age of the stellar population. For observed clusters with 
stellar masses below the minimum mass, the ionising flux is 
better accounted for by a single star, rather than by a popula- 
tion of massive stars, so that the use of synthesis models is not 
appropriate in these cases. We have also evaluated the disper- 
sion in the distribution of the effective rates of ionising photons 
at the minimum mass, and found it to be around 60%. This 
means that the initial mass of observed stellar clusters must be 
larger than a minimum value in the range 3x10"* M© to 3 x 10^ 
Mq, depending on age and metallicity, in order for evolution- 
ary synthesis models that assume a completely sampled IMF to 
be applied. We have also found indications suggesting that the 
emission of [S 11] is much less affected by sampling effects than 
that of other ions. The He n lines are in turn the most affected 
by sampling, especially during the WR phase, and should not 
be used to constrain the evolutionary status of stellar clusters. 

We also studied the correlations between different effective 
rates of ionising photons, finding in some cases correlation co- 
efficients close to one. This result agrees with the observational 
finding that Hn regions over a vast range of scales, ranging 
from those ionised by single stars to those ionised by super- 
stellar clusters, are found in a relatively narrow band in some 
emission-line diagnostic diagrams. 

Tables containing the effective rate of ionising pho- 
tons, their corresponding N, the corresponding LLL, and 
the complete correlation matrix for H'^, He°, He^, N°, 
N+, N++, O", 0+, 0++, S", S+, S++, C°, Ne", Ne+, Ar^, 
Ar^, and Ar^^ ions are available at the www address 
http : //laef f ■ inta . es/users/mcs/SED/ They include 
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the results using the CMHK code with standard and high mass- 
loss rate"^. 
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Fig. 2. Evolution of the cross-section- weighted number of ionising photons Qa^(X') as a function of tim e for different metal licities 
and ions, using Starburst99 (stellar tracks with enhanced mass-loss rates, and model atmospheres bv'S mith"et IDl|2QQ2l)). The 
values are normalized to 1 M© for a cluster with a Salpeter IMF in the mass range 0.1 - 120 Mq. The five main panels corre- 
spond to the metallicities Z^O.OOl (top-left), Z=0.004 (top-right), Z=0.008 (middle), Z=0.020 (solar, bottom-left) and Z=0.040 
(bottom-right). Each main panel is further subdivided into four sub-panels: the top-left sub-panel corresponds to H", He'', and 
H^; the top-right sub-panel to N", N^, and N^^; the bottom-left sub-panel to O'', O^, and O^^; and the bottom-right sub-panel to 
S°, S+, andS++. 




Fig. 3. Evolution of the cross-section-weighted number of ionising photons Qa^iX') as a function of time for different metallicities 
and ions, using the CMHK code (stellar tracks with standard mass loss rates). Panel division and symbols like in Fig.|2| 
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Fig. 5. Ai"""[Qa,(X')] values as a function of time fo r different ions and metallicities using Starburst99 with stellar tracks with 
enhanced mass-loss rates, and atmosphere models bv lSmith et alJ ll2002l The values are normalised to 1 M© for a cluster with a 
Salpeter IMF in the mass range 0.1 - 120 Mq. Panel division and symbols like in Fig.|2| 
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Fig. 6. Same as in Fig.|5l but using the CMHK code and standard mass-loss rate tracks. Panel division and symbols like in Fig.|2| 




Fig. 7. Relative dispersion cr(QaJ/Qa,. as a function of time for different metallicities and ions. Panel division and symbols like 
in Fig. 121 
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Fig. 8. Linear correlation coefficient p(X', Y-') as a function of time for different metallicities and ions. The five main panels corre- 
spond to the metallicities Z=0.001 (top-left), Z=0.004 (top-right), Z=0.008 (middle), Z=0.020 (solar, bottom-left) and Z=0.040 
(bottom-right). Each main panel is further subdivided into four sub-panels: the top-left sub-panel corresponds to pQi^, He"), and 
P(H'', He^); the top-right sub-panel corresponds to the correlation coefficients of ions with similar ionising potential; the bottom- 
left sub-panel corresponds to correlation coefficients related to the j] parameter; and the bottom-right sub-panel to correlation 
coefficients related to the [N n] 6584/HQr vs. [O m] 5007/^S diagram. 



